//==============================================================================
// LSGAlignment.cc
// 
// Linear Space Global Alignmnet. 
// 
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "SequenceAlignment.hh"

long limitingCaseCount = 0;


//------------------------------------------------------
// * GetLSpaceGAlignment
//------------------------------------------------------
// Recovers the best globally aligned subsequence of two 
// sequences.  
//
void
LSGAlignment::GetAlignment()
{
  // Allocate two columns for score table and two columns 
  // for intersect forwarding table ("c" table). 
  mScoreTable.Allocate(4,mRows); 
  GetAlignment(0,0,mRows-1,mCols-1);
}

//------------------------------------------------------
// * LSpaceGAlignment
//------------------------------------------------------
// Divide and conquer approach to alignment.  Finds a 
// single point on the optimum alignment (u,v) and
// recursively splits the alignment problem around that
// point.  When the sub-problem reaches the scale of aligning
// a single pair of letters (or letters and gaps), the 
// alignment is inferred directly from the shape of the 
// sub-problem (2x1, 1x2, or 2x2, for gap-letter, letter-gap, 
// and letter-letter respectively).  The full alignment is 
// simply the concatenation of the alignments from the left
// and right half of each sub-problem.  Here, I construct 
// the alignment by right-end concatenation on the mAlignS, and
// mAlignT vectors, which works out because recursion is depth first 
// on the left half then right.  
//
void 
LSGAlignment::GetAlignment(int x1,int y1,int x2,int y2)
{
  int w = x2-x1;
  int h = y2-y1;
  int u,v;
  char ans;

  // Sanity check
  if ((w<0) || (h<0)){
	 cout<<"Invalid range.\n";
	 exit(1);
  }

  if ((w < 2) && (h < 2)){
	 HandleLimitingCase(x1,y1,x2,y2);
	}else{
	  u = ((x2-x1)/2)+x1;
	  v = FindInterceptWithBestPath(u,x1,y1,x2,y2);
	  
	  GetAlignment(x1,y1,u,v);
	  GetAlignment(u,v,x2,y2);
	}
}

//------------------------------------------------------
// * HandleLimitingCase
//------------------------------------------------------
// Using push_back only works if recursion solves left 
// hand side before right hand side.
//
void
LSGAlignment::HandleLimitingCase(int x1,int y1,int x2,int y2)
{
  int w = x2-x1;
  int h = y2-y1;

  // Little verification that code is still running.  On Macs, this 
  // also has the effect of giving the finder some time for the mouse,etc.
  //limitingCaseCount++;
  //if (!(limitingCaseCount%1000)) cout<<"LC:"<<limitingCaseCount<<endl;
  
  if (((w==1) && (h == 1)) ||
      ((w==0) && (h == 0))){
	 mAlignS.push_back(mS[x2-1]);
	 mAlignT.push_back(mT[y2-1]);
  }else{
	 if ((w==1) && (h == 0)){
		mAlignS.push_back(mS[x2-1]);
		mAlignT.push_back('-');
	 }else{
		if ((w== 0) && (h == 1)){
		  mAlignS.push_back('-');
		  mAlignT.push_back(mT[y2-1]);
		}else{
		  cout<<"Invalid limiting case detected."<<endl;
		  cout<<"w="<<w<<" h="<<h<<"\n";
		  exit(1);
		}
	 }
  }
}

//------------------------------------------------------
// * FindInterceptWithBestPath
//------------------------------------------------------
// Finds the intersection of the optimum path with column
// u.  This uses the "c" table method described in Durbin, Eddy, 
// Krough, and Mitchison, "Biological Sequence Analysis".
//
int
LSGAlignment::FindInterceptWithBestPath(int u,int x1,int y1,int x2,int y2)
{
	int diagScore,leftScore,topScore,temp;
	int rows = y2-y1+1;
	int cols = x2-x1+1;
	int current = 1,previous = 0;   // score table columns
	int Ccurrent = 2,Cprevious = 3; // "c" table columns
	int i,j;

	// Adjust u to local sub-table coordinates 
	u = u - x1;

	// Fill in initial column with gap scores. 
	for(i = 0;i<rows;i++){
		mScoreTable[previous][i] = i*mGapScore;
	}

	// Compute optimum alignment score one column at a time for 
	// the given range.  
	int x,y,yIdx,xIdx = 0;
	do{
	  xIdx++;

	  // Fill in the top of each column with gap scores.
	  mScoreTable[current][0] = mGapScore*xIdx;

	  // Top row of "c" table can only come from top,left corner.
	  mScoreTable[Ccurrent][0] = 0;

	  for(yIdx = 1;yIdx < rows;yIdx++){
		 // -1 adjusts for offset between string and table coordinates
		 int  matchScore = Score((char)mS[x1+xIdx-1],mT[y1+yIdx-1]);

		 diagScore = mScoreTable[previous][yIdx-1] + matchScore;
		 topScore = mScoreTable[current][yIdx-1]  + mGapScore;
		 leftScore =  mScoreTable[previous][yIdx]  + mGapScore;	
		 temp=max(diagScore,leftScore);
		 temp=max(temp,topScore);
		 mScoreTable[current][yIdx]= temp;

		 // Update the "c" table which will map back the intercept.
		 // The "c" table is here implemented simply as two more 
		 // columns, columns 2 and 3, of the score table.  
		 if (xIdx == u){
			mScoreTable[Ccurrent][yIdx] = yIdx;
		 }else{
			if (xIdx > u){
			  int c;
			  // Note that the order of these if statements imposes 
			  // an order on alignment recovery.  The order given 
			  // prefers the diagonal (matches) first.  
			  if (temp == topScore) c=mScoreTable[Ccurrent][yIdx-1];
			  if (temp == leftScore) c=mScoreTable[Cprevious][yIdx];
			  if (temp == diagScore) c=mScoreTable[Cprevious][yIdx-1];
			  mScoreTable[Ccurrent][yIdx]=c;
			}
		 }
	  }
	  
	  Swap(previous,current);
	  Swap(Cprevious,Ccurrent);

	  if (!(xIdx%1000)) cout<<xIdx<<endl;

	}while((xIdx + x1) < x2);  	
	//}while(xIdx < (cols-1));  KJD ???

	// ScoreTable contains the traceback position for the -local- 
	// sub-table.  Add y1 to map it back to the full table. 
	// NOTE: I return Cprevious because the last swap made 
	// Cprevious the final column. 
	return(y1+mScoreTable[Cprevious][rows-1]);
}

